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model and the general properties of models displaying self-organized criticality. 
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I. INTRODUCTION 

Sandpile models were first introduced by Bak, Tang, 
and Wiesenfeld as explicit models of self-organized criti- 
cality (SOC) [3. Since then a vast literature has analyzed 
sandpile properties resulting from various definitions of 
the toppling rules; see Refs.^i2c4 for recent reviews. 

In this paper we present the relaxation and station- 
ary properties of a gradient-driven sandpile model with 
a metastable toppling rule. Interest in this type of 
work originates from seismology, where quasiperiodic be- 
havior of seismic activity has been observed for certain 
faults and investigated with SOC-related models The 
model we present can describe characteristics of systems 
undergoing rheologic flow, e.g., in the mining environ- 
ment or tectonic plates. In such systems stress can accu- 
mulate in various parts of the system for long periods of 
time, in contrast to normal fluid flow in which the local 
relaxation time is much smaller than the hydrodynamic 
time scale. 

Furthermore, we believe that the properties we de- 
scribe are of wider interest to other fields of nonequi- 
librium statistical mechanics where, e.g., stripe pat- 
terns similar to those that we observe appear in a va- 
riety of extended systems, including sand and biological 
systems I^Q- Our approach is quite general since it uti- 
lizes only a consistent local ordering of topplings, follow- 
ing the instantaneous maximum gradient, and the notion 
of a metastable site. Furthermore, the model we propose 
shows a quasiperiodic time behavior, a feature already 
mentioned for SOC-related models in "s] and explored 
recently in a similar context by Chapman 9] . 

The additional feature of the proposed model is the 
emergence of a spatial structure for the metastable con- 
figuration driven by the transport of grains through the 
system. A study of this spatial configuration in and close 
to the stationary state is the main objective of the fol- 
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lowing sections. 

The paper is organized as follows. In Sec. II we de- 
scribe the toppling rule and its connection with related 
models. In Sec. Ill we present the ID variant of the 
model and in Sec. IV the results for the 2D version, 
followed by conclusions in Sec. V. 

II. THE TOPPLING RULE 

We consider a gradient-driven sandpile with a toppling 
rule which takes into account not only whether a lo- 
cal threshold gradient is exceeded, but also whether this 
situation is the result of the addition of a grain to the 
site under evaluation. We introduce this rule as a sim- 
ple description of the dynamical weakening introduced in 
Ref.^:5]. 

The general description of the dynamics is as follows. 
Grains are dropped on the sites in a designated region 
of the lattice called the source zone; following a relax- 
ation rule, to be described in detail subsequently, grains 
can change their position on the lattice until they reach 
an open boundary and leave the system (lattice). In 
short, the model describes the transport of grains from 
the source zone to the open boundary. 

Let us now analyze in detail the toppling rule in 2D 
on a square lattice. At a given moment of time each 
site of the lattice has an associated height h of the grain 
column, and we also associate with it the set of gra- 
dients G = {gh9r,9u,gd}, where ga = h - and 
ha G {hi,hr,hu,hd} is the the height of the sand col- 
umn at the four nearest neighbor sites left, right, up, 
and down. 

We use two thresholds in our algorithm. (i)gmax is the 
stability threshold for the maximum of G. If at least one 
of the gradients of G is larger than gmax the site is called 
metastable. It topples only if this state is the result of the 
receipt of a grain, either from a neighboring site which 
had toppled, or when a grain was dropped on the site at 
the start of an updating run. Metastable states do not 
topple when a gradient larger than gmax develops as a re- 
sult of the loss of grains on neighbouring sites, (ii) gmin 
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is the minimum positive gradient, which fixes the condi- 
tion to stop the topphng; we call it the activity thresh- 
old. Once a site starts toppling, it sends grains, one at a 
time, along the instantaneous maximum gradient of G to 
its nearest neighbors, provided that max{G} > gmin- If 
there is more than one instantaneous maximum gradient 
a random choice is made among them. If one site topples 
it will send grains to some of its nearest neighbors. We 
refer to the process of relaxation of one site as a "top- 
pling" and to the sites that have received sand grains 
during toppling as "updated sites". Their coordinates 
are kept in a list for the next step of the dynamics. 

In one time step all the updated sites produced at the 
previous time step are checked for stability and, if im- 
stable, they relax according to the above algorithm. We 
specify that once a site is unstable and chosen to top- 
ple the algorithm will finish the toppling sequence at the 
site and then will move to another site. Physically, this 
is equivalent with neglecting the toppling time. 

Because we always topple along the instantaneous 
maximum gradient we can introduce a time order of the 
topplings. On physical grounds it seems reasonable to 
assume that the site that first receives a grain will topple 
first at the next time step, if it is unstable. We model this 
using a first in, first out list: the coordinates of the up- 
dated sites are stored sequentially in a list. In one time 
step the algorithm reads the list sequentially from the 
first entry: if the current site is metastable, it is toppled 
and the sites that are updated in this process are stored 
in the list for the next time step; if the current site is sta- 
ble, it is simply discarded from the list. We refer to this 
time ordered toppling rule by the acronym TOTR. (Re- 
call that a gradient toppling rule implies a non-Abelian 
sandpile where the order of toppling has definite conse- 
quences; see, e.g., Refs.0,0|.) 

Alternatively, we can disregard the time ordering and 
select at random a site from the list of updated sites 
for inspection of its stability. In Sec. Ill we make a 
comparative study of these two variants of the toppling 
rules in the 2D case. We refer to this randomly ordered 
toppling rule by the acronym ROTR. 

In other words, we can consider one time step of the 
dynamics, Af, divided into N bins {N being the number 
of sites of the lattice), each bin containing at most the 
coordinates of one updated site since in a time interval 
At/N we can assume that typically at most one site is 
active. An active site generates a set of updated sites 
when topples. In this description the TOTR fills the bins 
of the next time step in the order in which the updated 
sites are produced at the current time step; meanwhile 
the ROTR fills the bins of the next time step randomly. 
Evidently, most of the bins are empty as the number of 
updated sites is much smaller than the system size and 
the algorithm discards the empty bins using the list of 
updated sites. 

Thus, to resume, at a given instant of time we have 
three kinds of sites in our system: (i) inactive sites with 
all associated gradients less than gmax, (ii) metastable 



sites, which have at least one associated gradient larger 
than gmax] and (iii) active or unstable sites, which are 
the toppling sites at the given instant (they are perturbed 
metastable sites). 

One observation about the times scales in the model 
is now in order. We can think in terms of the existence 
of three time scales. The smallest one is the time in 
which a site has toppled, Tt- The second time scale is 
the surviving time of an updated unstable site, r^. The 
dynamics of our model is such that ^ tj. The third 
time scale is the time between two droppings of grains 
on the lattice Tq. After one grain is dropped the system 
relaxes globally through an avalanche of topplings. Here 
we are ensured that Ta 3> N^Ts, where Nl is the average 
number of time steps for the system to relax for a given 
size L. 

We define the size of an avalanche as the total number 
of topplings generated in a given relaxation process after 
one dropping. 

We stress that our toppling rule allows and intro- 
duces sites with unstable gradients after one avalanche 
has taken place. They emerge as the neighbors of the 
toppling (active) sites toward which the gradient is neg- 
ative. As an unstable site topples, the negative gradient 
increases in absolute value, but the sites along such a di- 
rection do not receive any grains and hence are not top- 
pled, according to our rule. When topplings in such an 
avalanche stop, such a site can accordingly have a max- 
imum gradient that is larger than the threshold value 
gmax- We call such a site metastable — it can sustain 
gradients larger than the threshold as long as it remains 
unperturbed. Physically, we associate this rule with a 
certain local metastability of the medium through which 
transport takes place. 

We make the observation that the metastable sites may 
be important in the regions where there is no dropping of 
grains, since they can live much longer than the average 
time between two droppings on the same site. 

As elaborated in the next two sections, this model has 
properties similar to the extended version of the forest fire 
model described in Ref.0, that is, quasiperiodic behavior 
in time and a spiked avalanche distribution. 

To link further with previous studies, it is also of in- 
terest to understand the relation between our proposed 
model and the condition for SOC in sandpile models. 
References |^ ^] present an analysis to establish the 
generic conditions a model has to satisfy in order to 
present SOC. The authors show that in two dimensions, 
or larger, and for conservative dynamics, an intrinsic spa- 
tial anisotropy is required to produce SOC for models 
which can be treated perturbatively. The model we pro- 
pose has a conservative toppling rule (the noise is also 
conservative in the region where no dropping takes place) 
and evolves in two spatial dimensions under conditions 
of anisotropy, but it does not present the features of SOC 
for the avalanche distribution. 

We think that the explanation of this anomaly resides 
in the fact that the stationary state of this model cannot 
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be characterised as a perturbation of an interaction-free 
model. As we present in full detail in Sec. IIVI the sta- 
tionary state of this model is characterized by the ap- 
pearance of deep and narrow valleys along the direction 
of grain transport. These features cannot be obtained as 
a perturbation of a diffusionlike relaxation. 

We also make the observation that in the proposed 
toppling rules there is no external noise term, e.g., sim- 
ilar to the noise term in the Langevin equation. The 
randomness in our model comes from the selection of 
the toppling order in the case of the ROTR, and from 
the random choice of the toppling direction in the case 
of a degenerate maximum gradient. Another source of 
randomness in these models is the 'internal noise', in the 
sense discussed in Ref. 0| , coming from the fluctuations 
of the many particle dynamics which may be present even 
in the case of pure deterministic microscopic equations. 
There is no straightforward connexion between the 'in- 
ternal noise' and the 'external noise' term of a Langevin 
equation. 

In connection with this theoretical aspect we mention 
that recently a model with metastable states ("sticky 
grains" is the term used by the authors) but with a 
height-driven toppling rule and dissipative dynamics was 
studied in Ref. [l^ and shown to be in the directed per- 
colation universality class. 

In the following sections we present a detailed discus- 
sion of the ID case and a statistical analysis of the men- 
tioned patterns of metastable states for the 2D case in 
the stationary state and close to the stationary state. 



III. THE ID CASE 

We start our study with a one-dimensional sandpile. 
We choose a lattice of dimension L with an open bound- 
ary aX X = L and with a wall at a; = 0. The grains are in- 
jected randomly in the region x G [1, w], with 1 < w < L, 
called the source zone. We choose the stability threshold 
dmax — 2 and the activity threshold gmin — 1- Let us 
analyze in detail the appearance of a metastable site in 
ID lattice. For simplicity we start with an initial stair 
configuration which has a monotonic step of maxjC} ~ 2 
descending from x = 1 to x = L. This configuration is 
marginally stable; if we drop a grain anywhere on the lat- 
tice an avalanche occurs. First we consider a restricted 
source zone where grains are only dropped at x = 1. A 
grain dropped at a; = 1 will then topple until x = L for 
the chosen configuration. At x — L an extra toppling 
takes place since gmin = 1, and the toppling of the dis- 
turbed state continues as long as the maximum gradient 
is larger than gmin- [Formally, the boundary condition is 
that h{L + l) = 0, h{Q) = oo.] When the avalanche stops, 
we therefore see that the gradient between the sites L — 1 
and L is 3, hence we have a metastable site at L — 1. Now 
if we drop a new grain at a; = 1, the same consideration 
shows that the metastable configuration moves to L — 2. 
The process continues until the metastable configuration 
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FIG. 1: The avalanche size distribution in ID at i = 128, 
256, 512, and 1024 with the fixed ratio w/L = 0.1. We 
notice the peak in the top right corner of the plot. 

reaches the site x — 1. After that a series of avalanches 
reconstructs the initial configuration, adding a grain at 
site L, then at L — 1, and so on. 

The characteristic time period of the system is con- 
trolled by the time in which the metastable site travels 
from the site L to the site 1. This kind of behavior is 
preserved if we use a source zone with w > 1 but with 
the ratio w/L small. The avalanche size distribution is 
characterized by a peak at the end of the distribution sup- 
port (see Fig. which results from the avalanches pro- 
duced while the metastable site is present in the transport 
zone between w + 1 and L. The smaller size avalanches 
are produced after the metastable site has reached the 
source zone. Figure ^ shows that for fixed ration w/L 
the avalanche distribution scales with L. 

We close this section with the observation that the 
distinction between the TOTR and ROTR is irrelevant 
in ID since typically there is only one updated site in the 
algorithm. 



IV. THE 2D CASE 

A. General features 

In this section we compare the behavior of the system 
in two dimensions with the simple and well understood 
behavior in ID. We choose for study a rectangular lattice 
of size L with an open boundary condition at x — L and 
a wall at a; = 1, with a periodic boundary condition along 
the y axis. A simulation starts with the initial condition 
specified by a uniform slope of size 1 per lattice step 
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FIG. 2: The patterns created by sites at the bottom of 
valleys for a 512 x 512 lattice during the transient regime 
(a),(d) and in the stationary regime (b), (c), (e), (f). On 
the left side the dynamics is TOTR and on the right side 
the dynamics is ROTR. The dots, often merged to lines, 
mark the sites which have at least two negative gradients 
larger in absolute value than pmax, i.e., at least two neigh- 
bors are metastable. (See the text for discussion of the 
symbols ©.) 

along the x axis, the height ai x — L being 0. The initial 
slope along the y axis is set to 0. We choose the stability 
threshold gmax = 4 and the activity threshold gmin = 1- 
(We have also performed simulations starting with an 
empty lattice; the stationary state is the same and only 
the transient regime is longer.) We use both toppling 
rules: (i) with the time ordered toppHng rule (TOTR) 
and (ii) with a random ordered toppHng rule (ROTR). 

The grains are dropped randomly at a; = 1, y S [1, L]. 
This source configuration allows us to explore the behav- 
ior of the transport region of the sandpile under minimal 
perturbation, since only one grain topples at the bound- 
ary of the transport region and the metastable sites do 
not have their average lifetime constrained by the aver- 
age time between two droppings on the same site; hence 
the structure of the transport zone is influenced only by 
a minimal flow of grains. 

The initial condition we start with can be viewed as 
a collection of interacting ID sandpiles oriented along 
the X axis, and consequently one might expect to see 
two-dimensional avalanches created by the interaction of 



FIG. 3: The structure of the bottom of the valleys as a 
function of the lattice size L. From top to bottom L = 
128, 256, 512, 1024; TOTR on the left side and ROTR on 
the right side. The source is on the left side and the open 
boundary on the right side of each panel. The vertical lines 
mark the open boundary. 



different metastable sites. 

Surprisingly we found that in the stationary state the 
sandpile develops a structure of valleys along the x direc- 
tion separated by terraces of sites in the metastable state. 
As one can see from Fig. |21 the valleys are not purely 
ID; instead they fluctuate slightly along the transverse 
direction and also show branched structures. Visual in- 
spection of Fig. 121 shows that the transverse fluctuations 
and branching arc more pronounced in the case of the 
TOTR algorithm. 

After the system reaches the stationary state we notice 
that in the case of the TOTR the valleys do not change in 
time except for small fluctuations which appear toward 
the open boundary. In the case of the ROTR the val- 
leys do change in time even in the stationary state. We 
illustrate this in Fig. |2 where the snapshots (c) and (f) 
follow after 10^ droppings on the snapshots (b) and (e) 
in the stationary regime. We observe that in the case of 
the TOTR the two configurations are almost identical, 
while for the ROTR there is a clear difference between 
the two configurations around the points denoted with 
the symbol 0. We return to this point when we analyze 
the cluster size distribution and the correlation length 
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FIG. 4; The profile of tlie sandpile along the source line 
X — 1 (continuous line) and along the line next to it, a; = 2 
(dashed line), for TOTR (left) and ROTR (right). We 
show one moment in time of the transient regime, (a) and 
(d) , while (b), (c) and (e), (f) are in the stationary regime. 
The coordinate is the height in number of grains. Note the 
change of the profile height with time, and the fact that 
the profile is unchanged in the stationary regime. 



behavior. 

Another observation which we can make from visual 
inspection of the stationary patterns, is that the struc- 
ture of the valleys changes with system size. From Fig. 
Owe see that for L — 128 the valleys are predominantly 
onedimensional. When the lattice size increases, trans- 
verse fluctuation and branches appear. 

We obtain further insight about the nature of the sta- 
tionary state if we analyze the time evolution of the 
height profile along z (perpendicular to the x — y plane) 
at the boundary of the source zone. In Fig. 0]we show 
the time evolution of the height profile along y for x = 1 
and X = 2. We see that in the transient state, for both 
toppling rules, the source zone (x = 1) has heights close 
to the heights of the x = 2 profile. In this configuration a 
grain dropped in the source zone will move more probably 
along the x axis, leaving the source zone in one or a few 
steps. In contrast, in the stationary regime we see that 
a completely different configuration arises. The average 
height of the source zone is significantly smaller than the 
average height of the x — 2 profile. A grain dropped 
in the source zone will therefore first travel along the y 
direction, until it meets a minimum which is connected 
with a valley, see Figs. El(b), and (e). In the stationary 
state the profile at the source remain unchanged for both 
the TOTR and ROTR as Figs. |31(b), (c), (e), and (f) 
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FIG. 5: The avalanche size distribution: (a) TOTR, 
(b) ROTR. The system size takes the values L = 
128, 256, 512, 1024. The plateau region and the first peak 
scale approximately with system size L. 



show. The data are taken from the same configuration 
presented in Fig. El 

We close the presentation of general features of the 
model with an analysis of the avalanche size distribution 
in the stationary state. FigureElshows for both dynamics 
that the size distribution is not power-law-like, but rather 
similar to the ID case (compare Fig. 0, although with a 
richer structure of peaks. We remark that a scaling pro- 
portional with lattice size L holds approximately for the 
plateau region and the first peak in the distribution, as 
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in the ID case. Since the valley length scales with L, we 
think that this feature of the avalanche size distribution 
is determined by the propagation of the avalanches along 
the valleys. 



B. Cluster characterization 

Having presented a visual description of the valleys in 
the previous subsection, we now consider a quantitative 
approach. For a numerical description of the patterns 
shown in Figs. El and |21 we concentrate on the clusters 
formed by the bottom of the valleys. We define a site to 
be at the bottom of a valley if it has associated with it two 
or more negative gradients larger in absolute value than 
the threshold, i.e., at least two neighbors are metastable. 

A set of bottom of the valley sites is said to form a 
cluster if they can be spanned by a path stepping only 
to nearest neighbors along the x ot y directions. 

We analyze the following quantities: the total number 
of clusters Nc{t) as a function of time, the total cluster 
mass Mcit) as a function of time, the longitudinal and 
transverse correlation lengths and ^j., and, finally, the 
cluster size distribution at a given time, N(s,t). 

We start with the time evolution of the total number of 
clusters Nc{t) and the total mass Mc{t) — J2s ^^i^'^)- 
We remind the reader that in this section a time step 
corresponds to a grain drop and the subsequent relax- 
ation, and therefore we neglect the relaxation time per 
avalanche. 

The data were collected over a single run starting with 
the initial condition specified at the beginning of this sec- 
tion. We have collected the data for one point at time 
intervals 5 x 10^, 10^, 10^, and 5 x 10^ for the system sizes 
L 128, 256, 512, 1024, respectively To eliminate the 
short time fluctuation we averaged at each point over a 
window of size 10'* steps from which we selected 100 mo- 
ments equally spaced. We mention here that we are con- 
strained to use single runs since the characteristic time 
in which the system reaches the stationary state scales 
approximately with the third power of the lattice linear 
size. 

Figure El shows that the total mass of the clusters and 
the number of clusters have a jump before the stationary 
state is reached. The jump is more pronounced in the 
case of the TOTR and is clearer for large system size. 
We observe that the time to reach the stationary state 
scales with the third power of system size. This fact can 
be explained if we assume that the average slope of the 
system converges to a nonzero value as the system size 
is increased; then the volume of the accumulated grains 
scales like L^, which determines the minimal number of 
time steps necessary to drive the system to the stationary 
state. 

In studies of this kind of model it is customary to test 
the data for scaling behaviour f{t, L) = L"' f [t / U') . Our 
plot shows that one can define a dynamical exponent 
z = 3 for the relaxation time, but a simple proportional- 




FIG. 6: Number of clusters (top) and their total mass 
(bottom) as function of time for the TOTR and ROTR. 
In ascending order the plots correspond to the lattice sizes 
L = 128, 256, 512, 1024. 



ity with the system size at a given power is not sustain- 
able, although for L = 512, 1024 the assumption holds 
acceptably (the corresponding lines are parallel within a 
good approximation). However, the scaling hypothesis 
does not hold for the system size L = 128, suggesting 
that boundary effects have a characteristic length of this 
order of magnitude. From Fig. Owe see that close to the 
the source zone the valleys do not present branching. We 
note also that close to the open boundary valleys do not 
form, which makes it plausible that the scaling is strongly 
affected by finite size effects. 

For a geometric charaterization of the clusters we de- 
fine the following correlation lengths, following Ref. [l3l |: 
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FIG. 7: Correlation lengths (top) and ^± (bottom) as 
function of time [ L=128 (+), 256(x), 512(*), 1024(0)]. 
Left panels are for TOTR right panels for ROTR. 



where ^ is the mean square distance along the longitu- 
dinal, X, or transverse, y, direction for the sites belonging 
to one cluster averaged over all clusters of size s, which 
can be expressed in the formula 

1 1 

Tig S 

Here the index i =||,-L indicates the longitudinal and 

transverse directions, Xn\ denotes the coordinate of the 
sites belonging to the ath cluster of size s, and is 
the number of clusters with size s. In other words, for 
each cluster from a given configuration of the lattice we 
compute its average square displacement from the center 
of mass in transverse and longitudinal directions with 
Eq. Next we average over all the clusters from a 

configuration with weights specified by Eq. (Q) . 

Figure [7| shows that the longitudinal correlation 
lengths jump before the stationary regime for system 
sizes L = 512 and 1024, similar to the behavior of the 
cluster number and the total mass( see Fig. O. 



FIG. 8: Histogram of cluster sizes at four moments of time. 
We take ti and t2 in the transient regime with Cl|(^i) < 
(,\\{t2) and ts, ti in the stationary regime with £,\\{ti) < 
^11(44). We mention that data are averaged around each 
time moment ti, i £ [1, 4], as described at the beginning of 
Sec. IV B. 

It is intriguing that the longitudinal correlation length 
shows chaotic behavior in the transient state as the sys- 
tem size increases, especially in the case of the ROTR, 
even though we have averaged out the short time fluc- 
tuations. On the other hand, the time evolution of the 
total number of clusters does not present fluctuations of 
the same relative magnitude (see Fig. . We conjecture 
that the large fluctuation of the correlation lengths in the 
case of the ROTR is related to a high frequency of coa- 
lescence and breaking up of clusters before the stationary 
state is reached. 

The transverse correlation presents the same kind of 
chaotic behavior. In any event, we see that > 1 (which 
is physically significant) only for the TOTR with L > 
512. Visual inspection of Figs. |21 and O shows that it is 
this situation in which a significant number of branchings 
arc observable. 

To elucidate this point further, we study, for both 
types of dynamics TOTR and ROTR, the cluster size 
distribution at various moments of time for the lattice 
size L = 1024. In Fig. |H1 we plot the cluster dis- 
tribution at two consecutive moments in the transient 
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regime {ti = 8.7 x 10^, t2 = 8.75 x 10^) corresponding 
to a minimum and a maximum of ^|| observed in Fig. 
13 that is, C||(^i) < C||(^2), and in the stationary state 
(ts = 1.84 X 10^, ti = 1.845 X 10^) with the same condi- 
tion Ciife) <i\\it4) ■ 

The plots show that in the transient regime the states 
with larger ^|| have a significant population of large clus- 
ters for both dynamics. In the stationary state in the case 
of the TOTR the cluster distribution changes very little, 
while for the ROTR the the large size tail presents ob- 
servable fluctuations. This observation seems to point to 
a signature for a stationary state of the ROTR in which 
large scale fluctuations appear a the cluster of valleys, 
as shown in Fig. [3 (We have also checked this fact for 
other configurations.) 

To summarize the two-dimensional behavior, we have 
studied the geometric and relaxation properties of the 
transport region of the sandpile. We found that as the 
system approaches the stationary state, clusters of deep 
valleys appear. Specifically, we explored the time evolu- 
tion for the total number of clusters, the total mass of the 
clusters, the correlation lengths (longitudinal and trans- 
verse) of the clusters, and the cluster size distribution at 
various moments of time. The main distinction between 
the TOTR and ROTR is the fact that the configurations 
obtained from the ROTR are strongly affected by fluc- 
tuations of the longitudinal correlation length, which we 
have shown to be associated with fluctuation in the tail 
of the cluster size distribution. 

The characteristic time to reach the stationary state 
scales with the system size as L^, with z w 3, but the 
magnitude of the observed quantities does not follow a 
clear power law type of scaling. The range of data we 
have does not allow us to decide if this is due to strong 
correction from boundary effects, or if scaling is genuinely 
broken. The setting in of a chaotic regime at large L 
complicates the analysis even further. 

V. CONCLUSIONS 

We have analyzed a gradient-driven sandpile model 
with local metastable states using two variants of the 



toppling rule: one which keeps the time order of the list 
of updated sites and one which selects a site from the list 
at random. We found that in 2D metastable sites gen- 
erate a rough landscape with deep valleys along which 
grains are transported. 

The valleys are organized in clusters with a pronounced 
development along the direction of grain flow (^|| ^ ^j_). 
For the time ordered toppling rule, we found that trans- 
verse correlations develop as the system size increases. 
Also, for large system size we observed chaotic behavior 
for the correlation lengths, associated with fluctuations 
in the size distribution of the clusters. 

The avalanches produced in the resulting valleys have 
one-dimensional properties and show non-SOC behavior. 
We believe that this is an example of a model in the 
strong coupling regime which escapes the classification 
made in Ref. 10] for models that can be treated pertur- 
batively. Yet we do not have an exact mapping between 
our discrete model and the continuum Langevin equa- 
tions used in Ref. ulS]. This may be another source of 
the discrepancy. At the moment this is a purely numer- 
ical study of a descriptive nature. To clarify the previ- 
ous uncertainties further analytical insight seems to be 
required. It would, e.g., be interesting to deduce a sim- 
ple mean field kind of equation capable of explaining the 
irregular time behavior at large L for the quantities de- 
scribed in the paper. 
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